
/* This program builds up a matrix starting with initial conditions, errors, and coefficients.
// matrices are indexed with column major indicies (go down column, then to next column, etc.)

//When calling from matlab, the syntax is:

// [X] = buildVAR(M, B, E, IC)

// X is the KxT output matrix
// M is the mu vector (Kx1)
// B is the transition matrix (KxK)
// E is the error matrix (KxT)
// IC is the initial condition (Kx1) */

/********************************************************************************************/

#include <math.h>
#include "mex.h"

/* Define Input Arguments */

#define Xout_  plhs[0]          /*KxT matrix*/
#define Min_   prhs[0]          /*Kx1 matrix (mu vector)*/
#define Bin_   prhs[1]          /*KxK matrix (transformation matrix)*/
#define Ein_   prhs[2]          /*KxT matrix (error matrix)*/
#define ICin_  prhs[3]          /*Kx1 matrix (initial conditions in)*/

void mexFunction ( int nlhs, mxArray *plhs[],
                   int nrhs, const mxArray*prhs[] )
                   
{ 
    double *IC, *M, *B, *E;
    double *X_ptr;
    int K, T;
    int t, i, k;
    double temp;
    
    T = mxGetN(Ein_);       /*get value of T*/
    K = mxGetM(Ein_);
    
    /*assign pointers to inputs*/
    M = mxGetPr(Min_);
    B = mxGetPr(Bin_);
    E = mxGetPr(Ein_);
    IC = mxGetPr(ICin_);
    
    /*set up output matrix*/
    Xout_ = mxCreateDoubleMatrix(K,T,mxREAL);
    X_ptr = mxGetPr(Xout_);
    
    
    /*built first period from initial conditions*/
    for (i=0; i<K; ++i) {
            temp = 0;
            for (k=0; k<K; ++k) {
                temp = temp + B[k*K+i]*IC[k];    /*temp is B*X0*/
            }
            X_ptr[i] = M[i] + temp + E[i];
    }
    
    /*now loop over the rest of the periods  (this loop starts at t=2)*/
    for (t=1;t<T;++t) {
        temp = 0;
        for (i=0;i<K;++i) {
            temp = 0;
            for (k=0; k<K; ++k) {
                temp = temp + B[k*K+i]*X_ptr[(t-1)*K+k];
            }
            X_ptr[(t)*K+i] = M[i] + temp + E[(t)*K+i];
        }
    }
                    
}

  
